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ABSTRACT 



^ ; We discuss the detectability of high-redshift galaxies via [Cn] 158/im line 

r v * ' emission by coupling an analytic model with cosmological Smoothed Particle Hy- 

drodynamics (SPH) simulations that are based on the concordance A cold dark 
matter (CDM) model. Our analytic model describes a multiphase interstellar 

l/-) 1 medium (ISM) irradiated by the far ultra-violet (FUV) radiation from local star- 

forming regions, and it calculates thermal and ionization equilibrium between 

t— i ■ cooling and heating. The model allows us to predict the mass fraction of a cold 

neutral medium (CNM) embedded in a warm neutral medium (WNM). Our cos- 
mological SPH simulations include a treatment of radiative cooling/heating, star 
formation, and feedback effects from supernovae and galactic winds. Using our 
method, we make predictions for the [Cn] luminosity from high-redshift galax- 
ies which can be directly compared with upcoming observations by the Atacama 
Large Millimeter Array (ALMA) and the Space Infrared Telescope for Cosmol- 
ogy and Astrophysics (SPICA). We find that the number density of high-redshift 
■ galaxies detectable by ALMA and SPICA via [C n] emission depends significantly 

on the amount of neutral gas which is highly uncertain. Our calculations suggest 
that, in a CDM universe, most [Cil] sources at z = 3 are faint objects with 
S u < 0.01 mJy. Lyman-break galaxies (LBGs) brighter than Rab — 23.5 mag 
are expected to have flux densities S u = 1 — 3 mJy depending on the strength 
of galactic wind feedback. The recommended observing strategy for ALMA and 
SPICA is to aim at very bright LBGs or star-forming DRG/BzK galaxies. 

Subject headings: cosmology: theory — stars: formation — galaxies: formation 
— galaxies: evolution — methods: numerical 
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1. Introduction 

Up to now, the high-redshift galaxies detected in large numbers are observed by the 
radiation from stars in rest-frame UV to optical wavelengths. Such studies have been fruitful 
in delineating the properties of star- forming galaxies at z > 3 that are bright in UV; i.e., the 
Lyman-break galaxies (LBGs) (e.g., Steidel et al. 1999; Adelberger & Steidel 2000; Shapley 
et al. 2001; Iwata et al. 2003; Steidel et al. 2003; Ouchi et al. 2004a,b). Within the hierarchical 
galaxy formation paradigm based on the cold dark matter (CDM) model (Blumenthal et al. 
1984; Davis et al. 1985), it has been shown that the spectro-photometric properties of LBGs 
can be accounted for if we associate them with relatively massive galaxies (M* > 1O 1O M ) 
situated in large dark matter halos (Mh a i > 1O 12 M ) (e.g. Mo & Fukugita 1996; Dave et al. 
1999; Nagamine 2002; Weinberg et al. 2002; Nagamine et al. 2004d, 2005b,a; Night et al. 
2006; Finlator et al. 2006). 

However, observing the UV light emitted by the young stars does not directly tell us 
about how much gas there is in the galaxy. In order to obtain a full picture of galaxy 
formation, we need to find answers to the questions such as: 1) How much neutral gas is 
available for star formation in high-redshift galaxies? 2) What is the baryonic mass (in 
particular, the neutral hydrogen mass) fraction as a function of halo mass, and how does it 
evolve as a function of redshift? 3) How does the volume-averaged neutral gas mass density 
evolve as a function of redshift? 4) What are the main cooling and heating processes of the 
ISM in high-redshift galaxies? 

One of the ways to detect the neutral hydrogen in high-redshift galaxies is to search for 
damped Lyman-a absorbers (DLAs, Wolfe et al. 2005, for a review) in quasar absorption 
lines. DLAs are defined as quasar absorption systems with a neutral hydrogen column density 
Ahi > 2 x 10 20 cm -2 , a threshold column density that effectively guarantees gas neutrality 
at high redshifts (Wolfe et al. 1986, 2005). Since DLAs are known to dominate the neutral 
hydrogen mass density at z > 3 (e.g., Lanzetta et al. 1995; Storrie-Lombardi & Wolfe 2000), 
it is expected that they represent a significant reservoir of neutral gas for star formation. A 
large number of DLAs have been discovered at z > 3, and they have proven to be one of the 
best probes of neutral gas in high-redshift galaxies, complementary to the optical to infra-red 
(IR) observations. A recent search for DLAs in the Sloan Digital Sky Survey (SDSS) data 
archive yielded a large sample of over 500 DLAs at z > 2.2 (Prochaska & Herbert-Fort 2004; 
Prochaska, Herbert-Fort, & Wolfe 2005). These observational results are beginning to put 
stringent constraints on theoretical/numerical models of galaxy formation. Furthermore, 
Wolfe et al. (2003b) opened a new window for probing star formation in high-redshift DLA 
galaxies utilizing the Cn* absorption lines, and this paper is motivated by their work. 

The Cn* absorption lines originate from the excited 2 Ps/2 state in the ground 2s 2 2p term 
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of C + . This is the same state that gives rise to [Cn] emission at A = 157.74 /xm through the 
2 Pz/2 P\j2 fine structure transition. Because this is the dominant coolant of diffuse ISM 
at temperatures T < 5000 K (Dalgarno k McCray 1972; Tielens & Hollenbach 1985; Wolfire 
et al. 1995; Lehner et al. 2004), detection of [Cn] emission is potentially another method for 
finding cold neutral gas at high redshifts (Petrosian et al. 1969; Loeb 1993). 

[C il] emission was first detected towards the inner regions of gas-rich spiral galaxies 
and starburst galaxies; i.e., dense star- forming gas irradiated by UV radiation from young 
star-forming regions near galactic centers (Russell et al. 1980; Crawford et al. 1985; Stacey 
et al. 1991; Carral et al. 1994). In this case, the brightness of the emission line suggests 
that it is produced in warm (T > 200 K), dense (% > 10 3 cm~ 3 ) photodissociated regions 
(PDRs) with pressures P/k B = 10 4 — 10 7 Kern" 3 at the interface between giant molecular 
clouds and fully ionized media. The [C II] luminosity accounts for 0.1 — 1% of the total far-IR 
luminosity of the nuclear regions. 

However, later space-based observations with the Long Wavelength Spectrometer (LWS) 
on-board ESA's Infrared Space Observatory (ISO) have shown that, for quiescent late-type 
galaxies like NGC6946, the [C n] emission from extended diffuse gas over the entire disk could 
also be significant compared to that from dense compact star-forming gas (e.g., Madden et al. 
1993; Malhotra et al. 1997; Leech et al. 1999; Malhotra et al. 2001; Contursi et al. 2002). 
The diffuse cold components of the ISM seems to be at least as important sources of the 
[Cn] emission as are compact regions (e.g., Sauty et al. 1998; Pierini et al. 1999; Pierini et al. 
2001; Contursi et al. 2002). 

More recently, [C n] emission has also been detected in the highest redshift SDSS quasar 
at z — 6.42 (Maiolino et al. 2005). Quasars are considered to host significant amounts of 
molecular gas (M gas pa lO n M ) and dust (Md ust ~ 10 8 M Q ), and are therefore forming stars 
at a significant rate of SFR > 1000 M yr" 1 (e.g., Omont et al. 1996, 2001, 2003; Carilli et al. 
2001; Walter et al. 2003; Beelen et al. 2004; Carilli et al. 2004). Given this large amount of gas 
contained in them, they should be very bright in [C n] luminosity, allowing them to be easily 
detected even at very high redshift. However, bright quasars have much lower comoving 
space density than LBGs and are not representative of the entire galaxy population; they 
are special sources hosting supermassive black holes, and may be in unusual dynamical states 
if their activity is triggered by mergers (e.g., Hopkins et al. 2005). 

The detection of [Cll] emission from normal (i.e., not quasars or active galactic nuclei 
[AGN]) high- redshift star-forming galaxies would be a more direct way than the absorption 
line technique to find out about the amount of neutral gas available for star formation. 
Such detections of [Cll] emission from high-redshift galaxies, for example, by ALMA and 
SPICA, would set an important constraint on the theory of galaxy formation. In particular 
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interferometric maps of [Cll] emission would reveal the spatial dimensions of the DLAs, 
a property that so far has eluded detection. By combining these measurements with the 
velocity widths of the [C n] emission lines, it will be possible to infer the dark matter masses 
of these objects, and the luminosity of the [C n] emission lines would tell us the heating rate 
of the gas. 

In this paper we study the detectability of high-redshift DLA galaxies by coupling an 
analytic model with cosmological SPH simulations that are based on the concordance ACDM 
model, and make predictions for upcoming observations by ALMA and SPICA. We also com- 
pare the computed [Cll] luminosity function with the observed LBG luminosity function, 
and discuss the connection between DLA galaxies and LBGs. Stark (1997) discussed the 
potential measurement of [Cll] emission from high-redshift galaxies, but he made various 
assumptions for the [Cll] luminosity function and the evolution of the characteristic lumi- 
nosity based on the local observations and numerical simulations, and focused more on the 
observing time and conditions rather than the intrinsic nature of high-redshift galaxies. The 
present study extends his original work, and is intended to be more physically motivated by 
using full cosmological hydrodynamic simulations which simulate structure formation from 
early epochs to the present time. We note that Suginohara et al. (1999) studied the de- 
tectability of atomic emission lines of carbon, nitrogen, and oxygen, but they focused on 
the intensity fluctuations owing to sources at z ~ 10 rather than the detection of individual 
objects. 



2. Cosmological Hydrodynamic Simulations 

The cosmological Smoothed Particle Hydrodynamics (SPH) simulations that we use in 
this paper were carried out using GADGET-2 (Springel 2005). This code adopts the novel 
'entropy-formulation' of SPH (Springel & Hernquist 2002), and includes 'standard' physical 
processes such as radiative cooling and heating by a uniform UV background of modified 
Haardt & Madau (1996) form (see Dave et al. 1999), star formation, supernova feedback, as 
well as a phenomenological model for feedback by galactic winds. The latter allows us to 
examine the effect of galactic outflows. (For now, we do not account for feedback associated 
with black hole growth, as in, e.g., Springel et al. (2005b, a); Di Matteo et al. (2005), but 
plan to investigate consequences of this mechanism in due course.) The ionization equi- 
libria of hydrogen and helium are solved assuming an optically thin limit throughout the 
simulation box. We will discuss this point in Section 4.2. We utilize a series of simula- 
tions of varying box-size and particle number to isolate the impact of numerical resolution 
on our results. The important parameters of the simulation runs are summarized in Ta- 
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ble 1. The same simulations were used for the study of the cosmic star formation history 
(Springel & Hernquist 2003b; Nagamine et al. 2004a), LBGs at z = 3 — 6 (Nagamine et al. 
2004d; Night et al. 2006; Finlator et al. 2006), damped Lyman-a systems (Nagamine et al. 
2004b, c, 2005c), massive galaxies at z = 2 (Nagamine et al. 2005b, a), and the intergalactic 
medium (Furlanetto et al. 2003, 2004, 2005). The adopted cosmological parameters of all 
runs are {tt M ,tt A ,Q b , a 8 , h) = (0.3,0.7,0.04,0.9,0.7), where h = H o /(100 km s" 1 Mpc" 1 ) is 
the Hubble constant. We also use the notation h 70 = h/0.7. 

The initial conditions for our simulations are set by tiny density fluctuations as mo- 
tivated by inflationary theories, and the calculations follow the development of structure 
using the laws of gravity and hydrodynamics. Therefore, in principle the cosmological simu- 
lations allow us to model galaxy evolution in a physical manner. However, even with current 
state-of-the-art computers, we still lack the computational resources to directly simulate the 
details of the ISM dynamics below ~ 1 kpc while at the same time solving for large-scale 
structure formation on tens of megaparsecs. Consequently, the simulations that we utilize 
in this paper do not directly resolve the multiphase nature of the ISM with temperature 
T < 10 4 K, but keep track of the dynamics of the hot ionized medium (HIM) with T > 10 4 K 
that is pressurized by the supernova feedback in star-forming regions, as we will describe 
later in Section 4.1. The mass of the cold gas with T < 10 4 K is then estimated with a 
sub- resolution ISM model as described in detail by Springel & Hernquist (2003a,b). In this 
sub-resolution ISM model, the hot and cold phase is assumed to be in pressure equilibrium 
(McKee & Ostriker 1977). The mass in high-density star- forming gas is dominated by the 
cold gas with T < 10 4 K, and the volume is dominated by the hot component. 

In the following sections, we will further supplement the simulation with a two-phase 
sub-resolution model in order to divide the cold gas with T < 10 4 K in the simulation into a 
cold neutral medium (CNM, with T ~ 80 K and n ~ 10 cm -3 ) and a warm neutral medium 
(WNM, with T ~ 8000 K and n ~ 0.1 cm -3 ). Therefore, in principle we are picturing a 
3-phase medium with HIM, WNM, and CNM (cf. McKee & Ostriker 1977). 

Within this multiphase ISM picture, Wolfire et al. (2003) estimated that the scale l P on 
which turbulent pressure begins to dominate the thermal pressure is ~ 215 pc for the WNM, 
based on the comparison of turbulent velocity dispersion and thermal velocity dispersion. 
Also the model of the vertical distribution of the Hi by Dickey & Lockman (1990) suggests 
that the FWHM of the vertical gas distribution of the WNM is ~ 500 pc (Wolfire et al. 2003). 
These two estimates indicate that the volume occupied by the warm ionized medium and 
the WNM could be substantial compared to that of the HIM in an ISM with a scale ~ 1 kpc. 
Furthermore recent high-resolution numerical simulations of a turbulent ISM suggest that 
the volume fraction of WNM could be comparable to that of HIM (Kritsuk & Norman 
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2002a,b; Kritsuk & Norman 2004). Since we are focusing on high-density ISM where the 
mass and density is dominated by the neutral component, here we restrict ourselves to the 
two-phase picture of the CNM and WNM, and in the following neglect the volume occupied 
by the HIM for the simplicity of the present sub-resolution model. Obviously this is an over- 
simplification of the actual 3-phase ISM, and could lead to an over-estimate of the volume 
fraction of CNM and thereby to an overestimate of the total [Cn] emission of galaxies. 
Therefore, the predicted amount of CNM gas based on our simplified model is more likely 
to be an overestimate rather than an underestimate. By the construction of our model, the 
total amount of CNM determined in our simulation cannot exceed the total amount of cold 
neutral gas (with T < 10 4 K, i.e., CNM+WNM) given by the simulation. 

3. The Model 

In this section, we describe a simple analytic model to estimate the mass fraction of 
the CNM for given mean gas density, FUV radiation intensity, metallicity, and dust-to-gas 
ratio. The notation that we use hereafter in the paper is summarized in Table 2. As noted 
earlier, our work is largely motivated by that of Wolfe et al. (2003b). These authors recently 
opened a new window for probing star formation in high-redshift galaxies utilizing the Cn* 
absorption line observed in DLAs. The idea is to use this absorption line to infer the cooling 
rate owing to the [Cn] fine structure line, which in turn tells us about the heating rate 
caused by the FUV radiation field from nearby star formation. One can then estimate the 
projected star formation rate (■?/>*) in DLAs as the heating rate is expected to be proportional 
to the local star formation rate (SFR). 

Wolfe et al. (2003b) considered a two-phase model in which the CNM and WNM are 
in pressure equilibrium. Subsequently Wolfe et al. (2003a) ruled out a solution in which 
Cn* absorption arises in the WNM, as it overpredicts the observed bolometric extragalactic 
background radiation. Howk et al. (2005) showed directly that Cn* absorption in one DLA 
cannot arise in the WNM. Since Cn* absorption arises from the same 2 Pz/2 state which gives 
rise to the [Cn] emission, we assume that the CNM is the dominant source for the [Cn] 
emission in this paper. 

3.1. Set-up & Key Equations 

Consider a ISM with mean density p and volume V in the vicinity of a star-forming 
region. This interstellar cloud consists of CNM and WNM that are in thermal equilibrium 
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within each cloud, and the two have density pc and pw, volume V c and Vw, respectively. We 
also denote the volume fraction occupied by CNM by fy = Vc/V , and the mass fraction of 
the CNM by fu- Given these definitions, the following simple relations immediately follow: 

PcVc + PwVw — PoV), (mass conservation) (1) 
Vc + Vw = V , (volume conservation) (2) 

, PcV c Pc , ,o\ 

Po Vo Po 

As mentioned above, in Equations (1) & (2), the volume of the HIM is ignored, and we are 
assuming that the volume fraction of the WNM is fairly substantial. Eliminating V w from 
Equations (1) & (2) gives 

fv = ^ = ALZ^l = 1 ~ (AW*) , and (4 ) 

Vq PC — PW ( £C_ _ I ) Pw_ 

\Pw J Po 

Po = (1 - fv)pw + fvpc- (5) 
Equation (5) can also be obtained from Equation (1) by dividing both sides by V . 
Combining equations (3), (4), & (5) gives 

, 1 1 ~ (pw/po) ((]) 

l - {pwlpc) ' 

Therefore, if we have estimates for pc/pw, and pw/po (or equivalently pc/pw, and fy), then 
we can compute the mass fraction of the CNM using Equation (6). Note that knowing the 
volume fraction of the CNM, fy, is equivalent to constraining the mean density po if Pc and 
pw are known, and the two quantities (fy and po) are related to each other via Equation (4) 
or (5). 



3.2. Phase Diagrams 

In order to identify the CNM and WNM densities, we compute the two-phase structure 
of a multiphase ISM by solving the equations of thermal and ionization equilibrium for a 
given gas density, FUV radiation field intensity, metallicity, and dust-to-gas ratio with the 
numerical techniques and iterative procedures outlined in Wolfire et al. (1995) and Wolfe 
et al. (2003b). We describe the further details of our method in Appendix A. 

Solving for the thermal balance results in a characteristic shape of pressure curves as 
a function of gas density as shown in the top panels of Figure 1 with a local minimum 



- 8- 



(P = P m in) and maximum (P = P max ) denned by <9(logP)/<9(logn) = 0. In this figure, 
the pressure curves are shown as a function of gas density for different input projected 
SFR per unit physical area (■?/>*) and metallicity (log(Z/Z Q )). Our full grid of models covers 
the range log^* [M yr" 1 kpc~ 2 ] = -4.0, -3.5, -3.0, -2.5, -2.0, -1.5, -1.0, -0.5, & 0.0, and 
log(Z/Z Q ) = —2.5, —2.0, —1.5, —1.0, —0.5, & 0.0. Two things can be read off from this figure: 
(1) When tp* — and only the UV background radiation is present, the resulting pressure 
curve appears almost identical to the case of logip* = —4.0; i.e., the impact of the UV 
background is negligible when log-?/'* 3> —4.0. (2) As ip+ (° r equivalently, the incident 
FUV radiation intensity) increases, the pressure maxima and minima move towards the 
upper right corner. For densities between those corresponding to P max and P m i n , the cooling 
curve is relatively insensitive to changes in temperature T. Therefore, at a given density, T 
undergoes a large increase when the cooling rate rises to match the increase in heating rate, 
thereby increasing the pressure. Conversely, outside this range in densities, the cooling curve 
increases rapidly with increasing T and as a result the pressure is insensitive to increases in 
heating rate. Also, as the metallicity decreases, metal-line cooling becomes inefficient, and 
the density and temperature have to be raised for a given heating rate, and therefore the 
equilibrium pressures go up. 

It is known that the two thermally stable phases can coexist over a narrow range of 
pressures P min < P < P max , where <9(logP)/<9(logn) > 0. Here, we follow the work by Wolfire 
et al. (1995) and Wolfe et al. (2003b), and assume that the equilibrium gas pressure equals 
the geometric mean pressure P = P geo = V PminPmax for identifying p c and pw, as indicated 
by the open triangles in the phase diagrams shown in the top panels of Figure 1. In passing, 
we note that Kritsuk & Norman (2002a,b) have shown, using high resolution numerical 
simulations, that the thermally unstable turbulent ISM settles down as a multiphase medium, 
where the majority of the masses (for both CNM and WNM) lie close to the local minimum 
pressure P m i n . 

For each equilibrium curve, the model of Wolfe et al. (2003b) gives the [C il] luminosity 
per hydrogen atom, £ c , as a function of gas density as we describe in Appendix A. The value of 
l c for the CNM can be determined at the CNM density as shown by the open triangles in the 
bottom panels of Figure 1. The value of l c increases with increasing metallicity as the heating 
of the gas (and correspondingly the cooling rate to match that) becomes more efficient via 
the grain photoelectric heating effect. At low metallicity (e.g., log(Z/Z ) = —2.5, the 
bottom left panel of Figure 1), the value of l c begins to decrease at high densities, because 
grain photoelectric heating becomes more inefficient owing to low dust content, and the 
heating becomes dominated by X-rays. The X-ray heating rate decreases with density since 
it becomes less efficient with increasing neutral gas fraction, which increases with density. In 
the case of log(Z/Z ) = —0.5, the value of t c continues to increase with increasing density as 
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the grain photoelectric heating remains efficient. The flattening of the curves at low density 
occurs when the cosmic microwave background (CMB) photons dominate the population 
rates of the fine-structure states. Note, the flattening is absent from the figure at very low 
metallicity since the population of the 2 Pz/2 state is proportional to metallicity when the 
CMB dominates. 



3.3. Mass fraction of CNM: f M 

Figures 2 and 3 illustrate the properties of our multiphase ISM model. Figure 2 shows 
PC-, Pwi the ratio of the two, and J'm as functions of projected SFR ip* for different metal- 
licities, while Figure 3 shows P m i n as a function of ip* for different metallicities. Figures 2a, 
2b, and 3 demonstrate that pw, Pc, and P m i n are smooth monotonically increasing functions 
of -0* as expected from Figure 1; as ^ increases, the equilibrium pressure P gco increases as 
does the equilibrium density. For a given tp*, as the metallicity decreases, the equilibrium 
density increases for cooling to match the heating rate. The CNM and WNM density scale 
roughly similarly, therefore the ratio pw/pc is relatively flat as a function of ■?/>* as shown in 
Figure 2c except for \og(Z/Z Q ) < —2.0 and log-?/** > —1.0. As a result, fu is a monotonically 
decreasing function of i/j*. Here, the values of /m were computed by Equation (6) assuming 
a constant value of fy = 7 x 10~ 4 , which is motivated by the following argument. 

Consider a column of length L and a cross section of unit physical area through the 
ISM. Let J'a denote the area covering fraction of the CNM, and n c \ the number density of 
the CNM cloud, and R the physical radius of each CNM (spherical) cloud. Then, 

An 

fv = —R 3 n ch and (7) 
f A = nR 2 n cl L. (8) 

Combining the above two equations and assuming typical values of /U ~ 0.5, L ~ 1 kpc, and 
R ~ 1 pc gives 

f v = 1.33 ~ ~ 1.33 x 0.5 x = 7 x 10~ 4 , and (9) 

L 1 kpc 

n c i = = 1.6 x KT 4 pc- 3 . (10) 

We take Ja ~ 0.5 because Wolfe et al. (2003b) find CNM gas in roughly half of the observed 
DLAs. If we instead adopt L = 100 pc, then f v = 7x 10^ 3 and n c \ = 1.6 x 10~ 3 pc~ 3 , which 
are in better agreement with the values adopted by McKee & Ostriker (1977). However the 
latter values result in greater than unity. 
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It is not surprising that the value of fy is very small, as the volume occupied by CNM 
clouds should be much smaller than that of the WNM. If we take Equation (9) and compute 
the values of fu using Equation (6), we obtain Figure 2d. The exact values of fy and fu are 
not important here, because when we later couple this model with cosmological simulations, 
we do not assume a constant value of fy. Instead, we take the mean density of the ambient 
ISM po directly from our hydrodynamic simulations and compute the value of fu for each 
gas element in the simulation using Equation (6). 

In summary, the procedure to obtain the CNM mass fraction in our model is as follows: 

1. Compute a grid of equilibrium curves in the phase space of pressure and density for a 
range of projected star formation rates and metallicities (top panels of Figure 1). 

2. For each equilibrium curve, identify pc and pw at P gco = V-Pmin-Pmax- 

3. Use Equation (6) and the inferred values of pc and pw to obtain fu- Here we take the 
mean gas density po directly from the result of cosmological simulations, and compute 
the amount of CNM gas for each dark matter halo as we show in Section 4.2. The 
model of Wolfe et al. (2003b) gives the [C il] luminosity per H atom £ c as a function 
of ip+ (see Appendix A), and we describe how we compute the [Cn] luminosity of each 
halo in Section 4.3. 



4. Results from Cosmological Simulations 

4.1. Star-forming gas in simulations 

Let us first look at the range of gas pressure and density that our cosmological hydro 
simulations cover. Figure 4 shows the distribution of cosmic gas in our cosmological sim- 
ulations in the plane of density vs. pressure. In the simulation, the gas that has density 
higher than the threshold density n t h = 0.13 cm -3 is treated as a multiphase medium and 
is able to form stars, as described by the sub-resolution ISM model of Springel & Hernquist 
(2003a,b). We consider this star- forming gas as the ambient ISM (with mean density po) that 
hosts CNM and WNM. The magenta dashed curve to the right of the SF threshold density 
is the analytic fit to the effective equation of state adopted in the simulation (Robertson 
et al. 2004). We stress that in Figure 4 we are plotting p Q vs. P for the simulations. The 
star-forming gas in the simulation is pressurized by heating owing to supernova feedback, 
and has a higher pressure than isothermal (P oc n) gas as traced by the magenta dashed line. 
The Q5, D5, and G5 runs also include a treatment of feedback owing to galactic winds. In 
this approach, a fraction of supernova feedback energy is given to star-forming gas particles 
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as kinetic energy and momentum in random directions, and some gas particles are ejected 
from star-forming regions owing to winds. Therefore, one can see that the spread of the gas 
distribution around the effective equation-of-state is larger in the Q5, D5, and G5 runs than 
in the 03 (no wind) run. The distribution of gas in the Q5, D5, and G5 runs is roughly the 
same. 

The three sets of two lines indicate the P m i n obtained from the pressure curves for 
metallicities of, from upper right to bottom left, \og(Z / Z & ) = —2.5 (black), —1.0 (blue), and 
0.0 (red), as a function of pc (the CNM density, shown in solid lines) and pw (the WNM 
density, shown in dot-dashed lines) at P geo , respectively. We stress again that in Figure 4 
we are plotting po vs. P for the simulations, whereas the analytic model calculation results 
show P vs. p c and p w . The criteria that we impose for the simulated gas to qualify for 
hosting CNM are P > P min and p > p w (see Equations [4] and [6]). The condition p c > pw 
is guaranteed by the construction of the model. 

The highest gas densities po attained in our simulations are n ~ 300 cm -3 in the Q5 
run, and n ~ 100 cm~ 3 in other runs. Most of the star-forming gas in the simulations has a 
higher pressure than P min in all runs, but we need to divide the gas into different metallicity 
ranges to perform a more detailed comparison. 

In Figure 5, the gas in the Q5 run is divided into different metallicity ranges, and is 
compared to the corresponding values of P min for the same metallicity range. It is seen that 
most of the star-forming gas, i.e., gas with n tn > 0.13 cm -3 , in the simulation is relatively 
metal-rich, and there is almost no star-forming gas with \og(Z/Z & ) < —2.5 that satisfies 
P > -Pmin- Since the simulation does not resolve the high density CNM directly, it is not 
a problem that the diluted density of the simulated gas does not fall in-between the range 
-Pmin < P < Pmax for the CNM densities pc- 

4.2. CNM gas in dark matter halos 

Before we compute the CNM mass in each dark matter halo, there are several steps in 
processing the simulation data. First, we identify dark matter halos by applying a conven- 
tional friends-of-friends algorithm to the dark matter particles in each simulation. After dark 
matter halos are identified, we set up a 3-dimensional cubic grid centered at the center of 
each dark matter halo covering the entire halo, with grid cell-size equal to the gravitational 
softening length of each simulation. Then, the gas mass, H I mass, metal mass, and star for- 
mation rate of each gas particle is smoothed over a spherical region of grid-cells, weighted by 
the SPH kernel. We now have all of the above quantities for each grid cell. These procedures 
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are the same as adopted by Nagamine et al. (2004b, c). 

We now compute the CNM mass contained in each dark matter halo. Given the four 
quantities (p , P,^, Z/Z Q ) in each grid cell that covers the dark matter halo and the two 
quantities (pc,Pw) obtained from the phase diagram shown in Figure 1, we can calculate 
the value of fu for each grid cell using Equations (6). Here, only those cells that satisfy 
P > Pmin an d p > p w (see Equations [4] and [6]) are considered to host CNM. Figure 6 
shows the volume-averaged value of Jm for each halo, (/m), where all the grid cells in each 
halo are equally averaged over and each point in this plot represents one halo. Contours 
are used to represent the scatter and avoid a saturation of points in the figure. The values 
°f (/m) are in the range 0.5 — 0.8 for relatively massive halos with Mh a io > lO n h~ 1 M Q . 
The 03 run has higher values of (f M ) owing to the absence of winds. The rapid decline 
of (/m) values at 10 < logM halo < 12 for the D5 and G5 runs owes to lower resolution in 
simulations with large box sizes. We consider that a true physical decline of (/m) should 
occur at logMh a i ~ 8.5, because at this mass-scale the DLA cross section decreases rapidly 
in the highest resolution runs as shown in Figures 2 & 3 of Nagamine et al. (2004c), and 
hence the CNM mass fraction is expected to decrease as well. 

Using the value of fu computed for each grid cell that covers the dark matter halo, we 
calculate the total CNM mass of each halo according to 

M C NM = Zif M p V , (11) 

where po and Vq are the gas density and the volume of each grid cell, and the index i runs 
through all grid cells for each halo. The resulting M C nm for each halo is shown in Figure 7. 
The shaded contours are for the CNM, and the black contour without the shade is the total 
neutral hydrogen mass in the halo. As we described in Section 2, ionization equilibria of 
hydrogen and helium are solved assuming an optically thin limit across the simulation box. 
However, in high density star-forming regions, the recombination time-scales are expected 
to be shorter than the ionization time-scales, rendering some reliability in our estimate of 
the total neutral hydrogen mass in each halo (see Section 5 for the prospects of improving 
upon this in the future). The CNM mass in our model cannot exceed the black contour by 
construction. The short- dashed, dash-dot, long-dashed lines indicate the 5%, 1%, and 0.5% 
of the halo mass, respectively. 

In the case of the 03 (no wind) run, the bulk of gas (both total neutral mass and the 
CNM mass) has mass fractions in-between 1 to 5% for most of the halos. These values 
are very close to that of the disk mass fraction 0.05 adopted by Mo et al. (1998). The most 
massive halos in the 03 run have slightly lower neutral mass fractions (~ 3%) than this value. 
When galactic wind feedback becomes strong (Q5, D5, and G5 runs), the total neutral mass 
fraction goes down to ~ 1%, and the CNM mass fraction to even lower values (0.1 — 1%). 
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This is because the gas is heated and ejected owing to winds, resulting in a lower neutral 
fraction. Similarly to Figure 6, the rapid decline of M CNM values at 10 < logM halo < 12 for 
D5 and G5 run owe to lower resolution in simulations with large box-sizes, and we consider 
that a true physical decline of M CNM should also occur at M ha i ~ 10 85 /i _1 M . 



4.3. [C n] emission from CNM 

Having obtained the CNM mass of each grid cell computed in the previous section, 
we can now compute the [Cn] luminosity of each halo by performing a similar sum to 
Equation (11) with an additional multiplication of [C n] luminosity per H atom, £ c , as follows: 

Lcn^UZ,^)^^, (12) 

where \i is the mean molecular weight and mn is the hydrogen mass. Notice that £ c is a 
function of metallicity Z and ip* of each grid cell. The resulting [Cn] luminosity of each 
halo is shown in Figure 8. Each point in this diagram represents a halo, but we use contours 
in order to avoid saturation in the scatter plot. The long-dashed and short- dashed lines 
correspond to the following scaling relationships: 

Lca = Cl (lW^Me) ' (13) 

where C\ = 10 41 and 10 40 ergs _1 , respectively. Most of the halos except for the least massive 
ones in the 'no-wind' (03) run follow C\ = 10 41 ergs -1 well, and the 'strong wind' (Q5, D5, 
G5) runs are better characterized by C\ = 10 40 ergs _1 . The decrease of L CII in the least 
massive halos in each run occur owing to same reasons as mentioned for Figure 6. 

Another way to characterize the distribution of [C II] luminosity of halos is to look at 
the cumulative luminosity function, which is shown in Figure 9. As expected from Figure 8, 
the 'no-wind' (03) run is brighter than the 'strong wind' (Q5, D5, and G5) runs by about 
an order of magnitude. The bright-end of the 03 and Q5 run is limited by cosmic variance 
owing to small simulation box-sizes, and the faint-end of the D5 and G5 runs are limited 
by the low resolution in large box-size simulations. The total luminosity function can be 
obtained by interpolating the results of the three runs when they overlap. We will perform 
such interpolation in the next section and in Figure 12. 
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4.4. Comparison with LBG Luminosity Function 



The connection between DLA galaxies and LBGs is of significant interest. The recent 
measurements of cross-correlation between DLAs and LBGs (Gawiser et al. 2001; Adelberger 
et al. 2003; Bouche & Lowenthal 2003, 2004; Cooke et al. 2006) suggest that the typical DLA 
halo mass could be similar to that of the LBGs' (M halo ~ 10 12 M ), and that DLAs could be 
strongly correlated with LBGs. Since LBGs contribute a significant fraction of star formation 
rate density at z ~ 3 and DLAs dominate the total Hi gas density at the same redshift, it is 
natural to expect that the two systems have some connection with each other. In this case, 
the energy input source for the heating discussed in Section 3 would be the central LBGs in 
the halo rather than the in situ star formation within DLAs (Wolfe & Chen 2006). 

In order to compare the computed [Cll] luminosity function with the observed LBG 
population, we show in Figure 9 the observed number density of LBGs at z = 3 with mag- 
nitudes Rab < 25.5 by Adelberger & Steidel (2000) and Adelberger et al. (2003) with a 
yellow shaded region. In addition, the following simple scaling laws and the observed cumu- 
lative luminosity function of LBGs can be used to obtain the magenta dot-long-dashed curve. 
Using the results of population synthesis calculations, Nagamine et al. (2005c) obtained a 
relationship between Rab magnitude and halo mass Afhaio &s 



where M halo is in units of h 1 M and C 3 = 55.03 (03 run) and 57.03 (Q5 run). Inserting 
this equation into Equation (13) gives 



For the case of the 03 run, this results in logL CII = —0AR A b + 51.01. Next, we com- 
pute the cumulative luminosity function of LBGs using a Schechter fit with parameters 
(m*,a,$*[/i 3 Mpc- 3 ]) = (24.54, -1.57, 4.4 x 10~ 3 ) obtained by Adelberger & Steidel (2000). 
We then convert the abscissa of this cumulative function using Equation (15) into [Cll] lu- 
minosity. As a result, the two magenta dot-long- dashed curves in Figure 9 are obtained: the 
brighter one on the right is for the 03 run scaling, and the fainter one on the left is for the 
Q5 run scaling. We also indicate the magnitudes Rab = 23.5, 25.5, 30.0, and 36.0 on this 
curve with filled triangles for the 03 run, which correspond to logLcn = 41.61, 40.81, 39.01, 
and 36.61, respectively. The agreement between the two magenta curves and the actual 
simulation results is not perfect, but the two magenta curves lie in-between the 'no- wind' 
(03) result and the 'strong wind' (Q5, D5, and G5) results, which we consider reasonable 
given the crudeness of the above scaling relationships. 



R AB = -2.51ogM halo + C 3 , 



(14) 



logL CII = -0.4 (Rab ~ C 3 ) + logd - 12. 



(15) 
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4.5. [C n] flux density 

The [C n] flux density of each halo can be computed by 

4ird L 2 4nd L 2 Au y ' 

where Au = v\^{y c jc)ri is the line width, z/i 58 = 1897GHz is the rest-frame frequency of 
the [C il] emission line, and v c is the circular velocity of the dark matter halo at a radius of 
overdensity 200, 



GM h ■ ^ 1/2 



'halo 



200 



1/3 /-, . x 1/2 



1/2 

(17) 



R 

- 123 Khw)''X^ km "' (18) 

The parameter rj relates the halo circular velocity to the velocity width of the line. In 
principle, there is a broad distribution of velocity width depending on the local dynamics 
and geometry of the line emitting gas, resulting in a large scatter of rj. Haehnelt et al. (1998) 
showed that, for DLAs, the median of the velocity width distribution is roughly 60% of the 
virial velocity of the halo. Here we adopt r\ = 0.6 following their work. Obviously this is an 
over-simplification of complex gas dynamics, but for now we are satisfied with this simple 
treatment. In the future we will compute the line profiles directly using the full velocity 
information in the simulations and study this issue in greater detail. 

Using the [Cn] luminosity shown in Figure 8 and Equation (16), we obtain the [Cn] 
flux density for each halo at z = 3 as shown in Figure 10. The distribution is roughly similar 
to that of [C n] luminosity vs. halo mass, and can be characterized by 

(M \ 2 / 3 

ira;) • < 19 > 

where C 2 = 10" 02 = 0.63 mJy (black long-dashed line, for the 03 run) and 10~ L2 = 
0.063 mJy (red short-dashed line, for the Q5, D5, and G5 runs). The flux density follows 
the scaling S v oc M^ Q , because Lcn oc Mh a i and Av oc v c oc M^J Q . We summarize these 
results in Figure 11. This figure shows the lowest limiting halo mass one can probe for a 
given flux density limit. The limiting halo mass is higher for the strong-feedback case than 
the no-feedback case, because the [C n] flux density is lower in the strong-feedback case for 
a given halo mass due to smaller amount of CNM mass. The dispersion around the scaling 
relationships were roughly estimated by eye to be ~ ±0.5 dex from Figure 10. 
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The cumulative flux density function is shown in Figure 12. Here, we combine the 
results of the Q5, D5, and G5 runs to cover the entire range of flux density as shown by the 
black solid line. The method for this interpolation is somewhat ad hoc, but the conclusion 
of this paper is not strongly affected by the details of this interpolation method. We simply 
require that the interpolated line (the 'Combined' result) smoothly connects the result of 
different runs. Similarly to the [C n] luminosity function, the runs with small box-sizes (Q5 
and 03 runs) are limited at the bright-end owing to a lack of massive systems, and the run 
with large box-size is limited at the faint-end owing to a lack of low-mass systems caused by 
the limited resolution. The result of the 03 run is brighter than that of the Q5 run roughly 
by an order of magnitude. Figure 12 suggests that the number density of galaxies at z = 3 
with flux density S„ > 0.1 mJy is about 2 x 10~ 2 Mpc~ 3 for the case of the 'no- wind' (03) 
run, and 3 x 10~ 4 Mpc~ 3 for the 'strong wind' run (the 'Combined' result). 

Like Figure 9, we show the observed number density of LBGs with magnitudes Rab < 
25.5 at z = 3 by Adelberger & Steidel (2000) and Adelberger et al. (2003) with the yellow 
shaded region. In addition, the two magenta dot-long- dashed curves are similarly obtained 
by using the same cumulative luminosity function of LBGs as in Figure 9. Inserting Equa- 
tion (14) into Equation (19) gives 

log S v = -0.27 {Rab - C 3 ) + log C 2 - 8. (20) 

For the case of the 03 run, this results in logS^ = —0.27 Rab + 6.7 (hereafter '03-scaling'). 
We use this scaling and obtain the two magenta dot-long- dashed curves in Figure 12: the 
brighter one on the right is for the '03-scaling', and the fainter one on the left is for the 
'Q5-scaling'. We also indicate the magnitudes Rab = 23.5, 25.5, 30.0, and 36.0 on this 
curve with filled triangles, which correspond to S u = 2.3, 0.7, 0.04, and 9.5 x 10~ 4 mJy, 
respectively, for the '03-scaling' case. The curve for the '03-scaling' roughly agrees with 
the actual simulation result, but the 'Q5-scaling' curve is in-between the actual 03 and the 
'Combined' results. 

Finally we look at the cumulative probability distribution function as a function of 
[Cll] flux density in Figure 13. The striking fact is that the fraction of the sources with 
S v > O.lmJy is very small at z — 3, i.e., less than 5%. The fraction of LBGs brighter than 
Rab = 30 mag is also about the same level. This is due to the large number of faint sources 
in the CDM universe, and it is simply a reflection of a steeply increasing number of dark 
matter halos with decreasing halo mass. Using the observed luminosity function basically 
gives the same result. The 03 run gives the most optimistic result in terms of the fraction 
of sources, which suggests that ~ 30% of all sources have S v > 0.01 mJy at z — 3. 
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5. Discussion & Conclusions 

We have coupled state-of-the-art cosmological SPH simulations with an analytic model 
of a multiphase ISM in order to compute the [C n] emission from galaxies at z — 3. We find 
that, in a ACDM universe, the majority of the sources are very faint with S u < 0.1 mJy This 
is presumably a generic prediction of the CDM model owing to an increasing number of low- 
mass dark matter halos with decreasing halo mass. If our model prediction on the faintness 
of the [Cn] sources is indeed correct, then it will be difficult for ALMA and SPICA to 
detect normal (i.e., not quasars, active AGNs, or strong submillimeter sources) high-redshift 
galaxies via [C II] emission in large number, as the sensitivity limit of both telescopes are 
expected to be S v ~ 0.1 — 1 mJy. The recommended observing strategy is therefore to focus 
on very bright LBGs with Rab ^$ 24 mag that are pre-selected through optical imaging and 
spectroscopically known redshifts. Since there is a large body of spectroscopic data on LBGs 
at z > 3, it should not be a problem to come up with observing targets. Our calculation 
shows that the brightest LBGs with R AB ~ 23.5 mag could have flux densities S v — 1 — 3 mJy 
depending on the strength of galactic wind feedback. 

Some caveats to note for the present work is that we have assumed that the physical 
conditions of neutral gas in DLAs and LBGs are the same, and that the dust-to-gas ratios of 
the gas are similar. If these assumptions are inappropriate for the LBGs at z ~ 3, then the 
above conclusions might not be entirely accurate. Our multiphase ISM model also assumed 
the SMC-type (silicates) dust composition which is less efficient for heating than other models 
such as the Galactic (carbonaceous) or the LMC-type models of dust grains. The true nature 
of dust in DLAs could be different from what we assumed here, and we expect a factor of two 
uncertainty in the dust-to-gas ratio. Note, however, that the metallicity of LBGs approaches 
solar, which is significantly higher than for DLAs. As a result, it is natural to expect that 
gas surrounding LBGs would have a dust-to-gas ratio higher than for the average DLA. 
Since the heating rate of the gas is proportional to dust-to-gas ratio, the [Cn] emission 
from LBGs is likely to be at the upper end of our calculation. In passing we note that 
star- forming DRG/BzK galaxies (van Dokkum et al. 2004; Daddi et al. 2004) would also be 
good candidates for [Cn] observation as they are considered to have higher star formation 
rates and dust content than LBGs and hence more luminous in [C n] luminosity. But their 
space density is lower by a factor of ~ 10 than that of LBGs. 

We have some confidence on the reliability of the properties of bright galaxies in the 
simulation, because we know from our previous work (Nagamine et al. 2004d; Night et al. 
2006; Finlator et al. 2006) that the simulations can reproduce the observed properties of 
LBGs at z = 3 — 6 relatively well. However the neutral gas is not well resolved near 
the resolution limit of each simulation, and this introduces some uncertainty for low-mass 
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galaxies (except for the Q5 run in which low-mass galaxies are well modeled). For this 
reason our scaling laws were mainly determined using the data for the massive galaxies 
in the simulation. There is also an intrinsic scatter in the CNM mass as a function of 
halo mass around the scaling relation (Figure 7). As a result, we expect 30% uncertainty 
in our predictions of the CNM mass fractions for high-mass halos, and somewhat larger 
uncertainties for lower mass halos. 

Even if the observation of normal high-redshift galaxies is difficult, the reward of [Cn] 
detection from such galaxies would be quite significant, because such measurements would 
directly constrain the amount of neutral gas and the properties of ambient ISM in high- 
redshift galaxies which is otherwise difficult to do. The observations of DLAs (e.g. Prochaska 
et al. 2005; Wolfe et al. 2005) have given us tremendous insights on the properties of neutral 
gas in high-redshift galaxies already, but [Cn] measurements will provide complementary 
information to further constrain theoretical/numerical models of galaxy formation. The 
combination of interferometric maps and the velocity widths of the [C il] emission could in 
principle tell us about the physical sizes of DLAs and the mass of hosting dark matter halos, 
as well as the physical properties of the gas such as the heating rate. 

This paper is just a first step towards such a goal, and a number of improvements in 
the analysis are needed in the future. For example, in the current simulations, an optically 
thin approximation was assumed throughout when solving for ionization equilibrium, and a 
simple radiative transfer calculation assuming a disk geometry was used to approximate the 
effect of radiative transfer from local star-forming regions. As computing power increases, a 
more accurate treatment of radiative transfer (e.g., Razoumov et al. 2005) and ISM physics 
on small-scales will become possible using direct information from simulations, such as star 
formation rates and metallicity of the gas. 

We acknowledge the significant contribution of Volker Springel to the simulations used 
in this work, and we thank for his useful comments on the manuscript. We are also grateful 
to Alexei Kritsuk for useful suggestions and discussions. This work was supported in part 
by NSF grants ACI 96-19019, AST 00-71019, AST 02-06299, and AST 03-07690, and NASA 
ATP grants NAG5-12140, NAG5-13292, and NAG5-13381. The simulations were performed 
at the Center for Parallel Astrophysical Computing at the Harvard-Smithsonian Center for 
Astrophysics. 
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A. Appendix: Details of the multi-phase ISM model 

Here, we briefly summarize the multi-phase ISM model developed by Wolfe et al. 
(2003b). The model is used for the calculation of Figure 1 in this paper. 

First, in order to relate the projected SFR per unit physical area, with the incident 
FUV radiation field J u , Wolfe et al. (2003b) assumed a plane parallel disk geometry with 
half-width h and radius R, in which the gas, dust, and stars are uniformly distributed. A 
disk is a reasonable approximation for dissipatively collapsing gas inside a dark matter halo, 
and it is also motivated by the observations of DLAs at z ~ 3 (e.g. Prochaska & Wolfe 
1997, 1998). In the cosmological simulations used in this paper, we still lack the numerical 
resolution (and/or necessary physics) to properly account for the formation of disk galaxies 
in correct number, which is known as the 'angular momentum problem' (Robertson et al. 
2004, and references therein). By assuming a disk geometry for the ISM model, we are 
implicitly assuming that the simulated galaxies also have the same geometry, which is a 
reasonable assumption. 



A radiative transfer calculation under this geometry gives 

+ 0(k u R) 2 ... (Al) 



Ju 2 Uvr 



when k u h <C k u R <C 1. The quantity is the luminosity per unit area on the sur- 
face of the disk, and k v is the absorption opacity of dust at frequency v. This equa- 
tion shows that J v depends on the aspect-ratio R/h only weakly. The relation Y, u = 
8.4 x 10" 16 (ip+ I M© yr" 1 kpc~ 2 ) [erg cm" 2 s _1 Hz" 1 ] is adopted from Madau & Pozzetti 
(2000). The effect of the UV background radiation is included in E„, but its effect is negli- 
gible compared to that owing to local star formation. 

Once the relationship between J v and ip* is obtained, we then compute the total heating 
rate T for a given value of ip* as follows: 

r = r d + r CR + rxR + r Cl , (A2) 

where the terms on the right- hand- side (RHS) are the heating rates owing to the grain 
photoelectric effect, cosmic rays, X-rays, and photoionization of C I by the FUV radiation 
field J v . In particular, the grain photoelectric heating rate per H atom, T^, is related to ■?/>* 
as 

T d oc ne J u [erg s _1 H _1 ] oc oc ip*, (A3) 

where k is the dust-to-gas ratio and e is the fraction of FUV radiation absorbed by grains 
and converted to gas heating (i.e., heating efficiency). The Yd was computed by adopting 
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the Weingartner & Draine (2001) expression for photoelectric efficiency in the case of pure 
silicates, blackbody FUV radiation, and extinction R v = 3.1. The heating rates T CR and 
T X R are also assumed to be proportional to ip*. 

The above heating rate has to be balanced with the total cooling rate A which includes 
the following terms: 

A = A FS + A MS + A LyQ + A GR , (A4) 

where the first term on the RHS is the cooling rate owing to fine-structure lines of [Cn] 
158 fxm (which dominates at T < 300 K) and [O i] 63 /zm (which is comparable to [C n] at 
T > 300 K), i.e., Afs = Ac n + Ao r The second term on the RHS owes to metastable 
excitations of C + , O , Si + , and S + , which becomes important at T > 3000 K. The third and 
fourth term are the cooling by Lja and grain recombination which become important at 
high temperatures. Cooling owing to transitions in the neutral species C°, Fe°, Mg°, and Si 
are not included as their contribution to A is negligible. 

We then let 

T = nA (A5) 

and solve the thermal and ionization equilibrium for given gas density n, dust-to-gas ratio, 
metallicity, and projected SFR. This leads to the two-phase gas structure as shown in the 
pressure curves in Figure 1. Once the thermal balance is achieved and CNM/WNM densities 
are identified, the [C il] emission per H atom can be obtained as 

4 = nA ClI +4 CMB , (A6) 

where the latter term on the RHS is the spontaneous energy emission rate in the limit of 
CMB excitation. The quantity £ c is shown in the bottom panels of Figure 1 as a function of 
density. 
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Table 1. Simulation Parameters 



Run 


Box- size 




WDM 


^gas 


e 


wind 


03 


10.00 


2 x 144 3 


2.42 x 10 7 


3.72 x 10 6 


2.78 


none 


P3 


10.00 


2 x 144 3 


2.42 x 10 7 


3.72 x 10 6 


2.78 


weak 


Q3 


10.00 


2 x 144 3 


2.42 x 10 7 


3.72 x 10 6 


2.78 


strong 


Q5 


10.00 


2 x 324 3 


2.12 x 10 6 


3.26 x 10 5 


1.23 


strong 


D5 


33.75 


2 x 324 3 


8.15 x 10 7 


1.26 x 10 7 


4.17 


strong 


G5 


100.0 


2 x 324 3 


2.12 x 10 9 


3.26 x 10 8 


8.00 


strong 



Note. - Simulations employed in this study. The box-size is 
given in units of Mpc, N p is the particle number of dark matter 
and gas (hence x 2), m DM and m gas are the masses of dark matter 
and gas particles in units of /i _1 M , respectively, e is the comoving 
gravitational softening length in units of h" 1 kpc. 



Table 2. Notation of Variables 



Variable 


Definition 


Pc, Pw 


density of CNM and WNM 


Po 


mean density of the total gas 


V C ,V W 


volume occupied by CNM and WNM 


V 


volume of the star-forming region under consideration 


/m 


mass fraction of CNM 


fv 


volume fraction of CNM 


!a 


area covering fraction of CNM clouds 


n c \ 


number density of CNM clouds 


R 


characteristic radius of the spherical CNM cloud 


L 


size of the star-forming region 



6 



CO 

I 



o c 
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X -25 
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o 
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III 1 
^log(Z/Z ) = -2.5 - 


i 1 i i i 1 i i i 1 i i i 1 

- log(Z/Z ) = -1.5 - 
MM 


i | i i i | i i i | i i i | 

- log(Z/Zj = -05 / - 
M l 


III 

l / f~7 /f~7 \ r~ 

- log(Z/Z ) = -2.5 - 

i 1 i i //^^^^^^^^^^^ 


- log(Z/Z ) = -1.5 - 






i 1 i i i 1 i i i 1 i i i 1 


j io g (z/z )=-o.5 ; 

i 1 i i M i i i 1 i i M 



-2 2 4 -2 2 4 -2 2 4 
log n [cm -3 ] log n [cm -3 ] log n [cm -3 ] 

Fig. 1. — Top panels: Phase diagram of pressure vs. density for different projected SFR of 
log^[M yr" 1 kpc~ 2 ] = -4.0, -3.5, -3.0, -2.5, -2.0, -1.5, -1.0.-0.5, and 0.0, from bottom 
(black) to top (red). The open triangles indicate the CNM and WNM densities at the 
geometric mean pressure P geo = V-fmax-fmin- Bottom panels: [C ii] luminosity per H atom, 
£ c , as computed by the model of Wolfe et al. (2003b, see Appendix A). See text for the 
qualitative trends of the curves shown in this figure. 
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Fig. 2. — Following quantities are shown as a function of projected SFR per unit physical area 
(ip* [M yr -1 kpc~ 2 ]): WNM density (log pw, panel [a]), CNM density (log pc, panel [b]), the 
ratio of the two (pw/Pc, panel [cj), and the mass fraction of CNM (/m, panel [dj) as computed 
by Equation (6) assuming a constant value of fy = 7 x 10~ 4 (see text). Different line types 
are for different metallicity as indicated in the legend of panel (c): log(Z/Z & ) = —2.5 (long- 
dash short-dash), —2.0 (long dashed-dot line), —1.5 (short dash-dot), —1.0 (long dashed), 
—0.5 (short dashed), and 0.0 (solid). 
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Fig. 3. — Projected SFR ^ vs. minimum pressure P min described in Section 3.2 and 3.3. 
Different line types are for different metallicities as indicated in the legend. The simulated 
gas is required to satisfy P > P min and p > Pw in order to host the CNM gas. 
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-6 -4 -2 2 4 -6 -4 -2 2 4 
log n [cm" ] log n [cm" ] 

Fig. 4. — Density vs. pressure of all cosmic gas in our 03, Q5, D5, and G5 simulations 
at z = 3. The vertical dashed line indicates the threshold density n t h = 0.13 cm -3 for star 
formation to take place in the simulation, and the gas to the right-ward of this threshold 
density is treated as ambient gas in star-forming regions in which the CNM and WNM are 
hosted. The 3 sets of 2 lines indicate the P min for metallicities of \og(Z/Z Q ) = —2.5 (black), 
— 1.0 (blue), and 0.0 (red) as functions of WNM (dot-dashed lines) and CNM (solid lines) 
density at P geo from upper right to bottom left, respectively. The lines for the WNM densities 
are shown because the criteria that we impose for the simulated gas to qualify for hosting 
CNM are P > P min and p > Pw ( see Equations [4] and [6]). The magenta dashed line is 
the analytic fit to the effective equation of state adopted in the simulation (Robertson et al. 
2004). The 4 contour levels are for (1, 10, 100, 1000) gas particles in each 2-dimensional bin 
of size (A logra, A log P/k) = (0.25, 0.29) from low to high level. 
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Fig. 5. — Density vs. pressure of the gas in the Q5 run at z = 3 for different metallicity 
ranges. The vertical dashed line is the same SF threshold density as shown in Figure 4. The 
two sets of lines indicate the P m ; n for the lower (blue) and higher (red) metallicity limit of 
each panel as functions of WNM (dot-dashed lines) and CNM (solid lines) density, similarly 
to those lines shown in Figure 4. The lines for the WNM densities are shown because the 
criteria that we impose for the simulated gas to qualify for hosting CNM are P > P m \ n and 
Po > Pw ( see Equations [4] and [6]). The majority of the gas in the simulation that satisfy 
the pressure criteria P > P min for the CNM is high metallicity gas. The contour levels are 
the same as in Figure 4. 
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8 9 10 11 12 13 9 10 11 12 13 
log M halo [h" 1 M@] log M halo [h" 1 M®1 

Fig. 6. — Mean mass fraction of CNM of each halo, (Jm), as a function of halo mass at 
z = 3. Each point on the figure represents the mean of Jm over all grid cells that cover 
each halo. The contours are used to avoid saturation in the scatter plot, and the four 
contour levels are for (1, 10, 100, 1000) data points in each two-dimensional bin of size 
(A log Mhaio, A (/a/)) = (0.1,0.02) from low to high level. The two highest contour levels 
are not seen well as there is a large pool of data points with (/m) = 0.0, particularly for 
low mass halos in each run. The rapid decline of (/m) values at 10 < logM ha i < 12 for the 
D5 and G5 run owes to lower resolution in simulations with large box-sizes, and we consider 
that the true decline occurs at around Mh a io ~ 10 85 /i _1 M Q (see text). 
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8 9 10 11 12 13 8 9 10 11 12 13 

log M halo [h 4 ivy log M halo [h- 1 ivy 

Fig. 7. — Shaded contours show the CNM mass of each halo as a function of dark matter 
halo mass at z = 3 in the simulations. The three contour levels are for (1, 10, 100) data 
points in each two dimensional bin of size (A logM halo , A logM CNM ) = (0.12,0.13) from low 
to high level. The black contour lines without the shade is the total neutral gas mass within 
each dark matter halo, i.e., the maximum amount of CNM that each halo could host. 
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Fig. 8. — [C n] luminosity of each dark matter halo as a function of halo mass at z — 3. Each 
point in this plot represents a single halo, but we use contours to avoid saturation of points 
in the scatter plot. The 3 contour levels are for (1, 10, 100) data points in each 2-dimensional 
bin of size (A logMhaio, A logMcNM) = (0.11,0.20) from low to high level. The long-dashed 
line in the top left panel and the short-dashed line in other panels show the relationship 
Lea. = G x (Afhaio/lO 12 ^- 1 ^©), where C x = 10 41 and lO^ergs" 1 , respectively. 
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Fig. 9. — Cumulative [Cn] luminosity function of simulations at z — 3. The ordinate is in 
units of comoving /iy Mpc~ 3 . The faint-end is truncated for the D5 and G5 run because 
of resolution limitations, and the bright-end of the 03 and Q5 runs is limited by cosmic 
variance owing to small simulation box-sizes. Note that the 03 run predicts a much higher 
[C n] luminosity than the Q5 run because of less efficient galactic wind feedback which allows 
more neutral gas to remain within the dark matter halos and emit [Cn] line radiation. 
The yellow shaded region indicates the observed number density of LBGs brighter than 
Rab = 25.5 at z = 3: n LBG = 4 x 1(T 3 h 3 Mpc" 3 (Adelberger & Steidel 2000; Adelberger 
et al. 2003). See text for the details on the two magenta dot-long- dashed curves, which are 
derived from simple scaling laws between halo mass and Rab magnitude of LBGs and their 
observed luminosity function. 
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Fig. 10. — [Cn] flux density of each dark matter halo as a function of halo mass at z — 
3. The 3 contour levels are for (1, 10, 100) data points in each 2-dimensional bin of size 
(A logM ha i , A log S„) = (0.11,0.13) from low to high. The long-dashed line in the top left 
panel and the short-dashed line in other panels show the relationship log S v = |(logM ha i — 
12) + C2, where Ci = —0.2 and —1.2, respectively. 
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Fig. 11. — This figure shows the lowest limiting halo mass one can probe for a given flux 
density limit, summarizing Figure 10. The dashed lines are the same scalings shown in 
Figure 10, and the shaded region shows the dispersion of ±0.5 dex at a given flux density 
around the scaling relation. 
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log [mJy] 

Fig. 12. — Cumulative [C n] flux density functions for 03, Q5, D5, G5 runs and the combined 
result of the latter 3 results. The ordinate is in units of comoving hj Mpc~ 3 . The rough 
detection limit of S v = 0.1 mJy for ALMA and SPICA is indicated by the vertical dotted 
line. The difference between the results of the 03 and Q5 runs owes to the difference in the 
strength of galactic wind feedback. The yellow shaded region indicates the observed number 
density of LBGs brighter than R AB = 25.5 at z = 3: n LBG = 4 x 10~ 3 h 3 Mpc~ 3 (Adelberger 
& Steidel 2000; Adelberger et al. 2003). See text for the details on the two magenta dot- 
long-dashed curves, which are derived from simple scaling laws between halo mass and Rab 
magnitude of LBGs and their observed luminosity function. 
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Fig. 13. — Cumulative probability distribution of [C n] sources as a function of flux density 
S v for the same models shown in Figure 12. It is seen that the majority of the sources are 
faint objects with S v < 0.1 mJy This suggests that one has to aim at very bright LBGs in 
order to have a detection even with ALMA and SPICA. 



